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Abstract 

Exact analytical results for the dynamics of two interacting qubits each of 
which is embedded in its own spin star bath are presented. The time evolution of 
the concurrence and the purity of the two-qubit system is investigated for finite 
and infinite numbers of environmental spins. The effect of qubit-qubit interac- 
tions on the steady state of the central system is investigated. 

1 Introduction 

Exactly solvable models play a very useful role in various fields of physics. They help 
improving our understanding of physical processes and allow us gain more insight into 
complicated phenomena that take place in nature Needless to recall the usefulness 
of exactly solvable models such as the harmonic oscillator, the nuclear shell model 
and the Ising model, to mention but a few. From a practical point of view, exactly 
solvable models serve as a very convenient tool for testing the accuracy of numerical 
algorithms, often used in the study of problems that cannot be analytically solved due 
to the complexity of the systems under investigation. 

In nature, quantum systems are influenced by their surrounding environment through, 
in general, complicated coupling interactions, leading them to lose their coherence [i2|. 
This refers to as the decoherence process IH [5]|. Moreover, quantum systems 
exhibit properties that do not have classical analogous 0. Of great interest is en- 
tanglement, the main ingredient for quantum teleportation and quantum computa- 
tion IHl |9l Uni im [El. Over the last years, many proposals have been made for 
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the implementation of quantum information processing. Solid state systems are very 
promising [fT3l [T4l and have been the subject of many investigations. In particular, 
decoherence and entanglement of qubits coupled to spin environments [15] attracted 
much attention fTG'.Tyj. Thus new exactly solvable models describing the dynamics of 
qubits in spin baths are highly welcome. Recently, the spin star configuration, initially 
proposed by Bose, has been extensively investigated [|T8l [T9i l20l |2T1 122]| . An exact 
treatment of the dynamics of two qubits coupled to common spin star bath via XY in- 
teractions is presented in [[23ll24ll . In this paper we propose to investigate analytically 
the dynamics when the two qubits interact with separate spin star baths. 

The paper is organized as follows. In section [21 the model Hamiltonian is intro- 
duced. In section [3] we present a detailed derivation of the time evolution operator 
and we investigate the dynamics of the qubits at finite N for some particular initial 
conditions. In section |4] we study the thermodynamic limit, in which the sizes of the 
spin environments become infinite. Section [5] is devoted to the second-order master 
equation. We end the paper with a short summary. 



2 Model 

The system under study consists of two two-level systems ( e.g., spin-| particles) each 
of which is embedded in its own spin star environment composed of N spins- 1. The 
central particles interact with each other through a Ising interaction; the corresponding 
coupling constant is equal to 45, where the factor 4 is introduced for later convenience. 
We shall assume that each qubit couples to its environment via Heisenberg XY in- 
teraction whose coupling constant is a, which is, in turn, scaled by A^^/^ in order to 
ensure good thermodynamic behavior. The spin baths will be denoted by Bi and i?2- 
The Hamiltonian for the composite system has the form 



where 
and 



H = Hq + Hs^Bi + HS2B2, (1) 
Ho = ASSlSl (2) 

TV TV 

Hs^B. = ^{SlY.S~' + S-Y.^+^^ (^ = 1,2). (3) 

k=l k=l 

Here and S'^ denote the spin operators corresponding to the central qubits, whereas 
S^'^ denotes the spin operator corresponding to the fc*^ particle within the i^^ environ- 
ment. Introducing the total spin operators J = J2k=i "^^^ ^^'^ ^ ~ XlfcLi ^'^^ of the 
environments B\ and B2, respectively, one can rewrite the full Hamiltonian as 

E = ASSlS''^ + -^{SlJ^ + S^_J+ + SlJ^ + Slj+). (4) 
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The dynamics of the two-qubit system is fully described by its density matrix 
p(t) obtained, as usual, by tracing the time-dependent total density matrix ptot{t), de- 
scribing the composite system, with respect to the environmental degrees of freedom, 
namely. 



p{t) = tTB^+B2[ptot{t)] 



U(t)ptot(0)Ut(t) 



(5) 



where U(t) and ptot(O) designate the time evolution operator and the initial total den- 
sity matrix, respectively. 

At t = the central qubits are assumed to be uncoupled with the environments; the 
latter are assumed to be at infinite temperature. This means that the initial total density 
density matrix can be written as 

Ptot(O) = p(0)®^® ^. (6) 

Here p(0) is the initial density matrix of the two-qubit system, and 1 is the unit matrix 
on the space C^®^. The former can be written as p(0) = J2ke Pke\Xk){Xi\^ with 

\Xe) £ {| ),\ — !")> I H — I + +)} for ^ = 1, 4. Similarly, we introduce the basis 

state vectors \ j, m) of C^®^, such that k < j < N/2 (k = for N even and k = 1/2 
for N odd), and —j<m<j. The time-dependent reduced density matrix can be 
expressed as 

p(t) = 2-2^5^p°,5^5^z/(iV,j>(iV,r)(j,r,m,s|U(t)|xfc)to|Ut(t)|j,r,m,s), 

k,l j,m r,s 

(7) 

where |j,r,m,s) = |j,m) ® \r,s), and iy{N,j) = (^/^J - (jv/2^j-i) Hence, 
our task reduces to finding the exact form of the matrix elements of the time evolution 
operator U(t) = exp{—iHt) (h = 1). This will be the subject of the next section. 



3 Derivation of the exact form of the time evolution op- 
erator 

The time evolution operator can be expanded as 

n=0 ^ ' n=0 ^ ' 

In order to derive analytical expressions for even and odd powers of the total Hamilto- 
nian H let us notice that i/o anticommutes with HsxB^ + Hs2B2^ that is, 

[Ho, Hs^Bi + Hs2B2]+ = 0. (9) 
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This can easily be shown using the following properties for spin-| operators: SzS± — 
±S±, and S±Sz = tS±. Moreover, it is easily seen that Hq"^ = 5^", which simply 
implies that for n > 0, 



(10) 



In the standard basis of ® C^, it can be shown that powers of Hs^Bi and HS2B2 
are given by 

[{J+J-f \ 

{J+J-f 

{J-J+Y 

\ {J-J+f) 

/ J+(J^J+)^ 



TT2k 



a \ 



2k 



(11) 



TT2k 

-"52B2 



^2k+l 
S2B2 



n) 













a \ 



2k 



AT/ 



J+(J_J. 


/ 

\ 





J.{J+J-f 

V 

({j+j-f 

{J-J+Y 

{J+J-f 

V {J-J+fj 

( 



000 J+{J-J+) 





V 



,(12) 

(13) 
\ 

a*) 



It follows that 

{HsiBi + Hs2B2Y^ 



E 

fe=0 



2/ 

^ 1 rriJB TT 
2^ )-"SlBl-"52i?2 



£-1 



2£ 



2fe „2{e-k) ^ ^ 
2£ 



E 

.A;=0 



fc=0 

£-1 



2£ 
2k + I 



Hft'HfB'^-' 



Dek + 



L 



tk 



(15) 



where 



Dik = diag 



(16) 



and 



Uk = an^^^i^a5[J+:^+(J_J+)'^(:^-J7+)'-'^-^ J+J'-(J-J+)^:r+J'_)'-*^-^ 
J_MJ^J_f{J_J^Y-^-\j_J^{J^J_f{J^J_Y-^-^ 



(17) 
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Using the fact that 



-1 



fc=0 

one obtains 



k=0 

2i 
2k + 1 



2e 

2k 



21 



2y/xy 



2i 



(18) 
(19) 



{Hsj^Bi + HS2B2 

I Ft 






\2£ 



a \ 



21 



X 



F2'- 



J+J- 



y/j-j+j+j- 



J-J+ 



where 



J-J- 



y/j+J-J+J- 



^2 



^3 



Ft 







Inserting equation (|20l) into equation (fTOl ). yields 

^2n 



21-1 



21 



2H 



2t 



1 

2 

( {M+Y + {M^ 




X 



J-J- 





(A4 + )"-(A^-)'^ 





{Mtr + {M,r 







J j jM+r-{M7^r 

{Mtr + (-M3 )" 









A 







(21) 
(22) 
(23) 
(24) 



(2)) 








(25) 



(Ml)" + (^4)" 
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where 







yj+j- ± VJ+J-) • 








Mf 






Mi 







The above operators satisfy 



{U^-^'U = (37) 



(26) 
(27) 
(28) 
(29) 



M± J+ = J+M3± , M± j; = J+M3±4, (30) 

M±J+J+ = J+J+M^, M^J+J^ = J+J-Mf. (31) 

Furthermore, one can show that the matrix elements of iJ^^+i ^re given by 

= ^6[{Mtr + (M^n (32) 

(i^^"+^)i2 = j^^-^^=[{Vj:j; + yTujimtr (33) 

+ (V:^- j+ - yj+jr)(M2 )"], (34) 



(^'"^')2i = J.^^^=[{^/j^+^/j^){Mtr (38) 

+ (v/:^-/4X)(A^rr], (39) 

{H'-+%, = -^'^[(A^^)" + (Al2)1, (40) 

(i^^"+^)23 = -(V2)J+J-^^^p^^^, (41) 

+ {./rji-^/j:j^){Mi)% (43) 
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2n+l> 



;3i 



-(5/2)J-j; 
1 



a/2 



VNJ-J+ 



(44) 
(45) 
(46) 

(47) 

(48) 
(49) 



2n+l^ 



Ml 



2n+l^ 



M3 



44 



(5/2) J-X 



a/2 



(7W+)" - (.Mr)" 

:[(v/:^+v/:7:x)(A<+)" 



2 ; j5 



J-- 



a/2 



(50) 

(51) 
(52) 
(53) 

(54) 
(55) 



Having in hand the explicit expressions of powers of the total Hamiltonian, it can 
easily be verified that the elements of the time evolution operator, obtained by inserting 
equations (|25l) and (|32l)- (|55l) into equation (O, are given by 



Uii{t) =^\cos{t\/ Mt ) + cos {t\/ )-i6 



sin (^t ^/Mf^ sin (^t ^/aT^ 



U2l{t) 



J- 



W2 rv/^+v/:^3m(t^ 



(56) 



, sin t\/A^i 

W2 r^+v/^3m(t^) 



(57) 



V J+J- - ^/ J+J- . ,, , . 

H — sm[t\ Ml 



(58) 
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iS 



1 



cos (^t \fMt^ ~cos{t yOwr ) 



]}. 



(59) 



(60) 



C/i2(i) 



ia/2 r^J+JL + v/X^ . 



(61) 



t/32W =^-^.2VJ.^- J-X {^"K^ V^) - ^Os(tVA^) 

sinf t 



+ i5 



2 ) sini J 



}. 



(62) 



ia/2 f x/J+JT + VXj; 



— J_ - < - — ^ , sm 



+ 



sm 



(63) 



1 r / / — 7\ / / \ r^^^(^v^^) sin(ty/A^ 

.-{cos(t^.M3") + cos(V.M3-)+^4 W 



C/iaW = - ^- 



sm(t/M+) 



C/23(i) 



+ i5 



1 



sm{ t\ M 



3 f' 



(64) 



(65) 



2,JJ_J+J+J. 
sinf t 



: I COS {t \j Mt^ - cos \J M2, j 



3 1 sinit Y •"''3 



(66) 



C/43(*) 



+ 



smU\/-M 



3 f' 



(67) 



f/24(t) 



COS (ty Mt^ + cos (ty j -i5 



sin [t^ Mt^ sm.{t^/M] 



^ za/2 fv/x;^+v/:o+ . / r— 

- J+ -.{- sin[t\ M 



sm{t\ M 



U: 



34 



^^^^-^^T^x^i — 7^1 — ^"i^ 



(68) 



(69) 



Ml 



Uuit) =J+J- 
-i5 



1 



(70) 



sml t\/M 



cos (^t \J -cos(t \Jm\ 



(71) 



It should be noted that the above components of the operator U(t) can also be 
derived by solving the Schrodinger equation fi22il 



For instance, we have 



df/n(t) _ 
dt 




dU2i{t) _ 


a 


dt 




dUsiit) _ 


a 


dt 




dUiiit) 


a 



a 



N 



a 

/AT 
a 



J.Un{t) - 6U2i{t) + -—J^U^,{t), 

V A* 



J_[/n(t) - SUsiit) + -^J^U.^it), 

V N 



J^U2l{t) + -^Xt/3l(t) +5f/4l(t). 

N V N 



dt 



(72) 

(73) 
(74) 
(75) 
(76) 



This set of differential equation can be solved by introducing the following transfor- 
mations: 



Unit) 



iSt 



Unit), 



U2iit)-^e-'''j^U2i{t) 



Usiit) 



iSt 



J-Usiit), 



f/4i(t) ^ e-'^'J-®J-Uiiit). 



(77) 
(78) 
(79) 
(80) 
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The resulting differential equations involve diagonal terms; they can be solved by tak- 
ing into account the initial conditions: 



f/.,(0) 



1 for i = j, 
for i ^ j. 



(81) 



Following the same procedure, it is possible to derive the remaining matrix elements 
of the time evolution operator. 

There exist many measures for entanglement. Here we shall use the concurrence, 
defined by ^ 



C{p) = max{0, 2 max[A/ Aj] — 



4 

E 

i=l 



(82) 



where the quantities Aj are the eigenvalues of the operator p{t) {ay ^ ay) p{t)*{(7y^ ay). 
The above measure is equal to one for maximally entangled states, and is equal to zero 
for separable states. The purity 



Pit) = iip{tf 



(83) 



can be used to quantify the decoherence of the central system; it is equal to \ for 
maximally mixed states, and one for pure states. 

It turns out that the density matrices corresponding to the initial product states 
|eie2), where ei = ±, are always diagonal. Furthermore, the numerical simulation 
shows that if the qubits are prepared in one of the above states, they remain unentangled 
regardless of the values of and 5. The purity decays less with the increase of 5. 



- — p(t) 

C(t) 



Figure 1 : The evolution in time of the concurrence (solid curve) and the purity (dashed 
curve) corresponding to the singlet state for 5 = a and N = 10. 



The matrix elements of the reduced density matrices corresponding to the states 

-^(|— +)±| +— )) and-^(| ++)±| )) are shown in the Appendix. The evolution in 

time of the concurrence and the purity corresponding to the above maximally entangled 
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Figure 2: The evolution in time of the concurrence (solid curve) and the purity (dashed 
curve) corresponding to the singlet state for S = Aa and = 10. 

states is practically the same; we only present the results obtained for the singlet state. 
It is found that, for fixed 6, the concurrence and the purity saturate as the number 
of spins increases. This naturally suggests the investigation of the case N oo 
(see the next section). For small values of the coupling constant 5, the concurrence 
decays from its initial maximum value Cmax = 1, then vanishes at a certain moment 
of time (i.e. entanglement sudden death [271). At long times, and sufficiently large 
N and 5, the purity and the concurrence converge to certain asymptotic values, which 
increase with the increase of the strength of interaction. Here it should be noted that, 
in contrast to the case of common spin bath, the singlet state is not decoherence free. 
This was expected because the latter state is not eigenvector of the Hamiltonian H. 
Nevertheless, we find that decoherence can be reduced with strong coupling between 
the qubits, in agreement with [22]. Finally let us remark that, although we only have 
considered infinite temperature, we can ensure that for long-range antiferromagnetic 
Heisenberg interactions within the baths, low temperatures will have the same effect 
on decoherence and entanglement of the qubits as strong coupling constants. 



4 Thermodynamic limit 

In the thermodynamic limit, ^ oo, the operators a/ J±J^/N converge to the posi- 
tive real random variable r whose probability density function is given by 

r ^ f{r) = Are-^''\ r > 0. (84) 



Indeed, it has been shown in |l22l |23l that the operator J+ / yN converges to the com- 
plex normal random variable z with the probability density function 



-e-^l^l^ (85) 
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Expressing z in terms of the polar coordinates r and 0, i.e., z = re^'^, simply gives 
\z\'^ = r^. Then integrating the corresponding probability density function over the 
variable (p from to 27r yields 

2 f 2 

dP{r) = f{r)dr = — l dcp r dre~'^^ 

J 



= Are-'^'' dr, (86) 

from which (|84l) follows. 

Hence we can ascertain that 

oo oo 

hm^2~^^tTB,+B,n(^V-J±-JT/N, VJ±Jt/n)= 16 j j rse-^^'''+''^n{r,s)drds, 



(87) 

where f2(r, s) is some complex- valued function for which the integrals in the right- 
hand side of equation (f87l) converge. 

Using the above result, one can express the nonzero elements of the reduced density 
matrix corresponding to the initial state -^(\ — h) — |H — the thermodynamic 
limit, as 

Pii(t) = p44(t) = ^[A+(t) + A-(t)], (88) 

P22(t) = p33(t) = ^[T+(t) + T_(t) + S+(t) + S_(t)], (89) 

P23{t) = -^[T+(t) + T_(t)+H+(t) + S_(t) + 2vI/(t)], (90) 
where ( we set a = 1 for the sake of shortness) 

oo oo 

(r±s? 



A±(t) = 16/ / rg e-'("'+^') , sin^it^S^ + {r ± sy)drds, (91) 

J J + {r ± s)^ \ / 



oo oo 

T±(t) = 16 f I rs e-2{^'+«') cos'{t^5'' + (r ± sf^drds, (92) 





oo oo 



S±(t) = 16y j rs e-^^^'+^'^ ^, ^ ^ sin^(tV^^ + (r ± g)2)drdg, (93) 



oo oo 

m{t) = j rs e-=^('^'+^')|cos(tV52 + + 5)2)003(^^^52 + _ 5)2 j 







^sin(^tv^52 + (r + s)2j sin(^t^/52' 



+ r - s 



+ ^ .2^/ ^ N2 V— 7 rfrrfs. (94) 

+ (r + sY O'^ + [r — sY 
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Unfortunately the above functions cannot be evaluated analytically; one should 
make recourse to numerical integration. This task can be significantly simplified by 
transforming the double integration into single one, which is much easier to carry out. 
To do that notice that the analysis of the expressions of the functions A±(t), T±{t), 
and leads to the evaluation of the probability density functions Q{fi) and R{ri) 

corresponding, respectively, to the random variables /i = r + s and rj = r ~ s (see ll28l 
for a similar situation). 

Let us begin with the variable fi; its probability density function is simply given by 
the convolution of /(r) with itself: 



Note that the upper limit of the integration over r is fi because the quantity fx — r should 
be positive. The evaluation of the integral is somewhat lengthy, but elementary; one 
finds that 

Qifi) = [2fi - v^e^'(l - 2/i2)erf(/x)]e-2'^', (96) 

where erf (x) designates the error function [|29l . 

Now consider the variable r] = r — s. One should be careful when using the 
definition of the convolution, since, in this case, t] belongs to the interval ] — oo, oo[. 
We have to distinguish between two cases, namely, ?7 > and r/ < 0. In the first case 
r G [0, cxd[, and hence 

Riv > 0) 



When T] <0, then r G [ 

Riv<0) = IQ j {7] + r)re~^'^''+'^^-^'^dr 

= l{-2r/ + v^e'''(l-2r/2)[l + erf(r7)]}e^2r,2_ ^^g^ 

Combining (l97l) and (l98l) . we obtain the following expression for the probability den- 
sity function of 7] over the real line: 

Riv) = l^vl + v^e^'ll - - erf(|r/|)]}e-2''^ (99) 

The above functions are depicted in figures [3] and IH Clearly, R{r]) is an even func- 
tion of its argument; it takes its maximum value at the origin, that is, raax{R{r])} = 



16 / (r7 + r)re-2('-+^)'-2'^'dr 







= i{2r/ + v^e^'(l-2r/2)[i_erf(r/)]}e-2'''. (97) 
-?7, oo[, which implies that 
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QO") 




-4-2 2 4 

Figure 4: The probability density function R{r]). 



R{0) = 0.886227. The maximum value of Q{fJ,) occurs at /iq = 1.14209, such that 
max{Q(/i)} = g(/io) = 0.859664. 

As a simple application let us prove the following: 
Theorem 1 The moments around origin of the random variables fi and t] are given 
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by: 



1 3 



l + 2"+^n2Fi(l + n,-;-;-l 



.2n+U r(f+l) 



2" 



^ + 2"(2n+l)2Fi(^ + n,i;|-l 



(r^2") = (/i2")-nv^r(i + n). 
= 0, 



(100) 

(101) 

(102) 
(103) 



where V[x), and 2-^1(0, 6; c; c?) denote the Gamma and the hypergeometric functions, 
respectively. 



Proof. Relation (|103l) is obvious since the function R{r]) is even. Let us prove (|100|) . 
We have that 



where 



2/n+l ~ -^n 



(104) 



(105) 



(106) 



To calculate Yn and /„, introduce the functions of the real variable x > 0: 



(107) 



IJx) 



00 



The first integral can be easily evaluated: 



(108) 







n+l 



(109) 



15 



The second integral satisfies 



Integrating by parts the RHS of (|108l) with respect to fi, and using (11091) . yield 

W(x) = ^^/.(x) + — 

Here we have used the fact that erf (x)' = 2e~^^/v/vr. 
Let = nlx'^'^^gnix). Then from (|1 1 1|) we have 

2(?i + l)^„,+i(x) = {2n + l)gn{x) + 



(110) 



(111) 



(112) 



On the other hand equation (IllOl) implies that 



X- 



dx 



+ {n + l)gn{x) = (n + l)5(„+i(x). 



(113) 



Combining the last two equations yields the following first order differential equation 
for the function gn{x): 



2x^i^ + gn{x) 



1 



0. 



dx " ^ ' + 
Differentiating both sides of (1114|) . and again using (|112|) . we obtain 



(114) 



d^ ^ /?, ^n + l\d ^ n + 1 
dx^ \2x x + l/dx 2a;(a; + 1) 



By setting y = —x, and hn{y) = gn{—x), we obtain 



d'^ / 3 n + l\ d n + 1 
+ — + t]^ + 



gn{x) = 0. 



(115) 



Idy'^ \2y y-lJdy 2y{y - 1 
which should be compared with the hypergeometric equation 



hn{y) = 0, 



(116) 



d'^ /c 1 + a + b — c\ 
dy'^ \y y-1 J dy ' y{y - 1) 



d ah 
+ 



iFi{a,h] c]y) = 0. 



(117) 



Thus 



n + 1, b 



2' 



3 
2- 



It follows that 



J„(x) =n!a;"+SFi(n + l,i;|;-x) 



(118) 
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Putting X = 1 yields 



771 

/„ = n!2Fi(n + l,i;|;-l), = (119) 



Also, using (lllll) , we obtain 



24+1 = (2n + l)nl 2Fi(n + 1, i; f ; -1) + (120) 



from which (1 1001) readily follows. The other moments can be evaluated with a similar 
method. ■ 

The functions (|9T1)-(|931) can easily be expressed in terms of the functions and 
R{r]). For example, we have: 

A+(t) = ^"Q(/i)^^^sin2(ty5^T7I^)rf/i, (121) 

A-(t) = r R{fi)^^^sm'(t,/6^TJ^)dfi. (122) 

It should be noted that in contrast to r and s, the random variables r] and /i are not 
independent. The function ^ (t) can not be further simplified, and should be evalu- 
ated using the double integration over the variables r and s. Nevertheless, using the 
Riemann-Lebesgue lemma, we can infer that 

lim ^{t) = ^{00) = 0. (123) 
In a similar way, the remaining functions tend asymptotically to: 

A+(oo) = - / Q{f^)^^^dfi, (124) 
^ Jo -\- fx 

A-(oo) = -j R{fi)^f^dfi, (125) 
T±(oo) = i (126) 
S+(oo) = - j Q{f,)^^—-^dfi, (127) 

S_(oo) = - j R{f,)^^—^dfi. (128) 



Notice that 



A±(oo) + S±(oo) = l (129) 
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independently of the values of 5. It follows that the asymptotic density matrix can be 
expressed as 

1 
Q 2-n 2-n 





p(oo) = 



4 

2-n 



2=& 

\0 









4/ 



(130) 



where 

It is easily seen that 



n = A+(oo) + A_(oo). 



limS+foo) = 0, limA-i-(oo) 

5^0 5^0 



1 

2' 



(131) 



(132) 



The corresponding asymptotic reduced density matrix reads 



n i 

p{oo)s=o 











1 

4 





o\ 






(133) 



which has a concurrence identically equal to zero. 

On the contrary, in the limit of strong coupling between the central qubits. 



lim S±(oo) = -, lim A±(oo) — 0. 

5— »oo 2 (5— »oo 



(134) 



Consequently, 



/O 





p{oo)s=oo = 

\0 

A straightforward calculation shows that 



1 

2 










(135) 



lim C(p(oo)) 

0— >oo 



1 

2' 



(136) 



In general, since < //^/ (/x^ + 5"^) < 1, then 



(137) 



This allows us to find the following explicit form of the asymptotic value of the con- 
currence: 

2-3n| 



C(oo) 



max 



{0. 



(138) 
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The latter can also be rewritten as: 



C(oo) 



2-3n 







for 
for 



2 
3' 

n < 1. 



< n < 

< 



(139) 



The variation of the asymptotic concurrence as a function of S is shown in figure [6l 
It can be seen that C (oo) remains zero up to a critical value Sc after which it increases, 
to tend asymptotically to |. The value of 6c can be evaluated numerically: 




Figure 5: The variation of C(oo) as a function of the coupling constant 5. The inset 
shows the critical point 5c. 



Sc = 0.342842, 



ni 



5=5c 



0.666667. 



(140) 



At the critical point, the density matrix reads 



Pc[00} 










1 

3_ 



1 
3 









(141) 
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5 Second-order master equation 

Under Born Approximation, the second-order master equation yields the following set 
of integro-differential equations: 

t 

hi{t) = -a' j (2pii(s) - p22(s) - P33(s)) cos[25{t - s)\ ds, (142) 



t 



t 

^13(0 = -a' J (2pi3(s)e^^^(*-) - P24(s)e^^^(*+^)) ds, (144) 


t 

p^^[t) = -a^ j 2pi3(s) cos[2(5(i - s)] ds, (145) 


t 

P22(t) = -O? j (2P22(S) - Pll(s) - P44(S)) COs[25(t - s)] ds, (146) 







t 



^23{t) = y" 2p23(s) cos[2(5(t - s)] ds, (147) 



t 

p24{t) = -a^ J (2p24(s)e2^^(-*) - pi3(s)e-2^^(*+^)) ds, (148) 



t 

h3it) = -o? j (2p33(s) - pii(s) - p44(s)) cos[25(t - s)] ds, (149) 



t 

^34(i) = / (2p34(s)e^^^(^-*) - Pi2(s)e-^^^(*+^)) ds, (150) 







Pu{t)^-a'^J (2pu{s) - p22{s) - p33{s)Jcos[2S{t- s)]ds. (151) 







Some of the above equations can be solved under a time-local approximation for 
which the matrix elements Pij{s) are replaced by Pij{t). One can find that (S and t 
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given in units of a and a respectively) 



pn(t) = ijl + [-1 + 2(p?i + pI,)] exp{i[cos(25t) - 1]} 

+ 2(p°, - p°J exp{^[cos(25t) - 1]} 
P22{t) = ^1 1 + [-1 + 2{pl^ + pOg)] exp{l[cos(25t) - 1]} 
+ 2(p^2 - P33) exp{^[cos(25t) - 

+ [-1 + 2(p°3 + pI,)] exp|l[cos(25t) - 1]} 
+ 2(p^3 - P22) exp|^[cos(25t) - 1]} 



(152) 



(153) 



P33(i) = T< 



(154) 



P44(t) = ^<ll + 



-1 + 2(p°4 + p?j] exp|i[cos(25t) - 1]} 



Pl4(t) 
P23(i) 



+ 2(p°4 - pJi) exp{^[cos(25t) - 
P?4exp{^[cos(25t)-l]}, 



P23exp<i ^[cos(2(5t) - r 



(155) 

(156) 
(157) 



These solutions describe approximately the dynamics at short times. In fact, the 
smaller the coupling constant 5, the better these solutions are. 

Note that when 5 = ( i.e. nonlocal dynamics), then 



exp{-^[cos(25t) - 1]} 



-2*2/" 



n 



1,2. 



(158) 



Thus the second order time-local master equation shows that the nonlocal dynamics, 
or, in general, the short time behavior follow a Gaussian decay law. Note that the solu- 
tions corresponding to the diagonal elements reproduce their asymptotic limit, namely, 
pii{oo) = |. However, those corresponding to the off-diagonal elements fail to repro- 
duce the steady state, since, for example, equation (|157l) implies that p23{t) 0. To 
end our discussion let us remark that equations (11431) . (11441) . (11481) and (11501) can be 
analytically solved only when 5 = 0. For instance (see figure |7]), 



Pi2(t) = ^[(p?2 + P^4)e-*^/^ + (p' 
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(159) 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 

at 

Figure 6: The variation in time of the the matrix element pn{t) corresponding to 
the singlet state. The solid curve represents the exact solution, and the dashed curve 
represents the approximate solution (I152I ). The parameters are iV = 10 and 6 = a. 

6 Summary 

In summary we have investigated the dynamics of two qubits coupled to separate spin 
star environment via Heisenberg XY interactions. We have derived the exact form of 
the time evolution operator and calculated the matrix elements of the reduced density 
operator. The analysis of the evolution in time of the concurrence and the purity shows 
that decoherence can be minimized by allowing the central qubits to strongly interact 
with each other. The short-time behavior, studied by deriving the second-order master 
equation, is found to be Gaussian. The next step may consist in considering more 
central qubits, and investigate whether the above results still hold. 

Appendix 

Using trace properties of the lowering and raising operators, it can be shown that the 
nonzero matrix elements corresponding to the initial maximally entangled states 4^(1 — 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 

at 

Figure 7: The variation in time of the the matrix element puit) corresponding to 
the singlet state. The solid curve represents the exact solution, and the dashed curve 
represents the approximate solution (I159I) . The parameters are = 10 and 5 = 0. 

+) ± I H — )) are explicitly given by: 

Pu{t) = 2~(^^+'hTB,+B,{Ul2ml2it) + Ul3mkt)}, (160) 

P22(t) = 2-^^''+'hTB,+B,{U22mUt) + U23mUt)}, (161) 

P23{t) = ±2-^'''+'hTB,+B,{U22mL{t)}, (162) 

P33(t) = 2-^^''+'hTB,+B,{Us2mUt)+Us3mUt)}, (163) 

Pu{t) = 2-(2^+^)trB,+B,{f/42(t)?7i2(t) + ?743(t)t/i3(t)}. (164) 
Those associated with the initial state ^(| ) ± | + +)) read: 

pn(t) = '2-^'''^'hTB,+B,{Uu{t)Ul{t) + U,,{t)Ul{t)y (165) 

P22(t) = 2-^^''+'hTB,+B,{U2imUt) + U2,{t)UUt)}, (166) 

pu{t) = ±2-^^''+'hTB,+B,{Un{t)Ul{t)y (167) 
P33(t) = 2-'-'''+\TB,+B,{u3i{t)Ul{t)+Uu{t)Ul{t)Y (168) 

P44(t) = 2-(2^+l)trB,+B.{f/4l(t)f/il(t)+f/44(t)f/i4W}- (169) 
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